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Abstract — We study the optimal control of propagation of 
packets in delay tolerant mobile ad-hoc networks. We consider a 
two-hop forwarding policy under which the expected number of 
nodes carrying copies of the packets obeys a linear dynamics. We 
exploit this property to formulate the problem in the framework 
of linear quadratic optimal control which allows us to obtain 
closed-form expressions for the optimal control and to study 
numerically the tradeoffs by varying various parameters that 
define the cost. 

Index Terms — Linear quadratic control, Delay Tolerant Net- 
works, Two-Hop Relay Routing. 

I. Introduction 

In DTN (Delay Tolerant Network) mobile ad-hoc networks, 
connectivity is not needed any more and packets can arrive at 
their destination thanks to the mobility of some subset of nodes 
that carry copies of a packet. A naive approach to forward 
a packet to the destination is by epidemic routing in which 
any mobile that has the packet keeps on relaying it to any 
other mobile that arrives within its transmission range. This 
leads to minimization of the delivery delay at a cost, however, 
of inefficient use of network resources (in terms of memory 
used in the relaying mobiles and in terms of the energy used 
for flooding the network). The need for more efficient use 
of network resources motivated the use of the more economic 
two-hop routing protocols in which the source transmits copies 
of its message to all mobiles it encounters, but where the 
latter relay the message only if they come in contact with 
the destination. Furthermore, timers have been proposed to 
be associated with messages when stored at relay mobiles, 
so that after some threshold (possibly random) the message 
is discarded. The performance of the two-hop forwarding 
protocol along with the effect of the timers have been evaluated 
in [1] which brings up the possibility of optimization of the 
choice of the average timer duration. 

The optimization of the two-hop relay routing as well as of 
extensions that we propose here are the central objectives of 
our present work. Instead of using fixed parameters, however 
(e.g. timers which have fixed average durations whose values 
can be determined as in [1]), we propose a dynamic optimiza- 
tion approach based on optimal control. Thus the parameters 
of the two-hop relay protocol associated with a message of a 
given age are allowed to change with time. 

DTNs attracted recently the attention of the networking 
community [5], [6], [7]. Among DTNs, a relevant case is 
that of mobile ad-hoc networks, including systems where 
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human mobility is used to diffuse information through portable 
devices [6]. 

Lack of persistent connectivity makes routing the central 
issue in DTNs [12], [17]. The problem is to deliver messages 
to destinations with high probability despite the fact that 
encounter patterns of mobile devices are unknown. Message 
diffusion algorithms trade off message delay for energy con- 
sumption, e.g., number of copies per delivered message and/or 
transmission range. As mentioned before, the naive solution 
is epidemic routing [18]. 

The two-hop routing protocol considered here was intro- 
duced by Grossglauser and Tse in [10]; the main goal there 
was to characterize the capacity of mobile ad-hoc networks 
and the two-hop protocol was meant to overcome severe 
limitations of static networks capacity obtained in [11]. Two- 
hop routing, in particular, provides a convenient compromise 
of energy versus delay compared to epidemic routing: the stan- 
dard reference work for the analysis of the two-hop relaying 
protocol is [9]. Fluid approximations and infection spreading 
models similar to those we use here are described extensively 
in [19]. Interestingly, we show (Remark 12. U that the fluid 
mean field approximation turns out to be an exact description 
of the dynamics of the expectation of the system's state due 
to the linearity of the dynamics under two-hop routing. 

Algorithms to control forwarding in DTNs have been pro- 
posed in the recent literature, e.g, [15], [8]. In [15], the 
authors describe an epidemic forwarding protocol based on 
the susceptible-infected- removed (SIR) model [19]. They show 
that it is possible to increase the message delivery probability 
by tuning the parameters of the underlying SIR model. In [8] 
a detailed general framework is proposed in order to capture 
the relative performances of different self-limiting strategies. 
Finally, under a fluid model approximation, the work in [2] 
provides a general framework for the optimal control of the 
broad class of monotone relay strategies, i.e., policies where 
the number of copies do not decrease over time. It is proved 
there that optimal forwarding policies are of threshold type. 

In this paper, we consider non-monotone relay strategies for 
two-hop routing, and apply optimal control theory to capture 
general trade-offs on message delay, energy and storage. After 
presenting the general model in the next section, we formulate 
the control problem in Section[In]and then provide the solution 
in Section [TV] In Section [V] we present and solve a similar 
control problem defined on an infinite horizon. An extension of 
the initial control problem is studied in Section[VT] Section[VTI] 
presents a numerical exploration, and Section IVIIII concludes 
the paper. 



II. The Model 

Consider K classes of mobiles, where class k has Nk 
mobiles. Let TV be a if dimensional column vector whose 
fc-th entry is N^. The time between contacts of any two nodes 
of respective classes i and j is assumed to be exponentially 
distributed with some parameter A, 3 Q The validity of such a 
model in the special case of a single class has been discussed 
in [9], and its accuracy has been shown for a number of mo- 
bility models (Random Walker, Random Direction, Random 
Waypoint). Let A be the K x K matrix whose ij-th entry is 
Xij . We assume that the message that is transmitted is relevant 
for some duration r. We do not make any assumption on as to 
whether the source or the mobiles know whether the messages 
have reached the destination within that period or not. 

Let the source be of class s and the destination of class d. 
Let £(t) be a ^-dimensional column vector whose j-th entry 
denotes the size of population of class j that has the packet. 
Let X(t) = E[£(t)]. We assume that each component of X 
evolves according to 



dXjjt) 
dt 



A si (N i -X i (t))-M i (t)X i (t),i = 1, 



,K, (1) 



where M,*(t) > 0. The above dynamics is a generalization 
of the well-known dynamics of two-hops routing (see for 
example [2]), containing the additive linear control term on 
Xi\ this new term represents the effect of timeouts by which 
message copies are discarded at intermediate relays of class 
i. We should point out that, as it will be clear later, these 
dynamics of different components of X will in fact be coupled 
by the introduction of a control term which will be picked 
optimally, as one minimizing a particular cost function that 
involves the entire vector X. In what follows, we now employ 
an equivalent expression for ([T} which simplifies the analytic 
derivation, i.e., 

dXi(t) 



dt 



= A si N-A id X i (t)-M i (t)X i (t),i = 1, . . . ,K (2) 



where Mi(t) = (A si — A id + Mj(t)). For future use, we 
introduce the compact notation 

dX(t) 



dt 



A m N - A out X{t) - MX[t) 



A, 



where we define Ai n — diag(A s i, . . . , A s k) 
diag(A ld , . . . , A Kd ), and M =diag(M 1) . . . , M K ). 

Assume that the source has a message at time and let T d 
denote the amount by which the message is delayed. Let "Jt be 
the CT-algebra generated by £ s , s < t. Denote the conditional 
successful delivery probability by := P(T d < t\J t ). 

Given 3^, the number of arrivals at the destination during 
time interval [0, t] is a Poisson random variable with parameter 
Y.^i^id fs =Q ti{ s )ds. Therefore 

K 

*(t) = 1 - exp ( - V A id I Zi(s)ds ) (3) 



E A - / ; 



'Note that our multidimensional description of the system allows one in 
particular to extend exponentially distributed inter contact times to the much 
larger class of phase type distributions. 



(An alternative more detailed derivation is given in the Ap- 
pendix.) Since exp(— x] is concave, we have by Jensen's 
inequality 



K 

E[*(t)] > D(t) := 1 - exp ( - ^ A td / 

i=i ^ s 



id i Xi(a)da) (4) 

s=0 



Remark 2.1: We note that related models using linear dif- 
ferential equations have frequently been used in DTNs to 
approximate the dynamics of the properly scaled number of 
mobiles that have a copy of the file; this is the limit of the mean 
field dynamics. This type of approximation has been shown 
to be tight in various related models as the total number of 
nodes N tends to infinity see e.g. [3]. As we just saw, it turns 
out that in the special case of two hop routing, the dynamics 
of the mean field limit coincides with that of the expectation. 
On the other hand, in the mean field limit, the inequality in 
© becomes an equality. 

III. Controlling the timers 

We may control the timers by allowing M to vary in time. 
Let u(t) = -M(t)£(t) and u(t) = E[u(t)]. We then have: 

dX{t) 



dt 



A m N -A out X(t)+u(t) 



where < X(t) < N, u(t) < , 



(5) 
(6) 



with the vector inequalities interpreted componentwise. 

The cost. The following performance measures could appear 

as part of the overall objective: 

1) Delay cost: We wish to maximize the lower bound given 
in <j4j on the success probability, i.e. on the probability 
that the delay T d does not exceed the threshold r beyond 
which the message is considered irrelevant. Hence we 
wish to maximize X(t), or some monotone function of 
it, where X(t) := £\ =1 A jd f* =Q Xj(s)ds. 

2) Indirect delay cost: We may wish to include a penalty 
for large values of u (which correspond to timers that 
expire at a high rate) as higher values of u may result 
in longer delays. 

3) Memory cost: We may wish to minimize or bound 
the number of copies in the system in order to avoid 
saturation of the memory available at the mobiles. This 
can be done by directly including a cost on each of 
the components Xj(r) where Xj(t) := j Q Xj(s)ds (or 
having an instantaneous penalty on Xj(t)). 

4) Energy cost for the network: The energy cost for class 
j mobiles is given by £ Aj d Xj(t). 

In view of the above, we introduce the following cost function 
corresponding to an initial state x and a policy u: 



V(t;x,u) = - Cl X(tf + c 2 



K 

E 



(uj(s) — u) 2 ds 



K K 



where a's are all positive. Let R := — ciAdAj + C3l+C4A^ ut , 
where is the column vector whose i-th entry is given by 
Aid- Then we can write 

V{t;x,u) = Vc 2 / (ui(s) - u) 2 ds + X(t) T RX(t) (7) 

We assume henceforth that i? is positive semi-definite, and we 
can also take C2 = 1 without any loss of generality. 

The objective is then to minimize V(t;Q,u). 

By state augmentation, we have a standard optimal control 
problem with quadratic cost (0 and linear state dynamics: 

^P- = -A out X(t) + u(t) + A m N, = X(t) 

at at 

The state and control are restricted according to (O. 

One could use the theory of constrained linear quadratic 
control, such as [4]. Or, one could choose to track a value 
that is sufficiently far from the boundary so that a controller 
without the constraints (as those in (O) will be satisfied in 
practice with a high probability. 

IV. Solution to the Optimal Control Problem 

Let Z (t ) : = (X, (t ) , . . . , X K (t ) , X x (t ) , . . . , X K (t ) ) T . Then, the 
composite state dynamics can be written as: 



dZ 
~dt 



Z = 



AZ + Bu + c, where A 



X 
X 



, B = 







-A out 

/ 

A m N ' 




and the cost function (expressed in terms of Z, and with 
terminal time r) becomes: 

J(t;Z,u) = Z t (t)QjZ(t) + / (u - u) T {u - u)dt 

Jo 



where with R > defined as earlier, Q f = 
Letting w := u — u, we can rewrite the above as 







R 



dZ , / A in N + u 

— = AZ + Bw + c where c := n 
dt \ U 



J(r; Z, w) = Z t (t)Q } Z{t) + / w T wdt 



Hence what we have is an affine-quadratic optimal control 
problem [14]. For each fixed finite r, this problem admits a 
unique strongly time-consistent optimal solution: 

w{t) = -B T (PZ + k) 

where Pit) > 0, kit) are unique continuously differentiable 
solutions of 

P + PA + A T P - PBB T P = 0, P(r) = Q f 
k + A T k + Pc-PBB T k = 0, fe(r) = 
and, min J = Z(0) T P(0)Z(0) + 2fc T (0)Z(0) + 2m(0) 



where m is the unique solution of 



1 



k T c~ -k T BB T k = 0, 
2 



m(r) = 



With this solution at hand, one can of course readily obtain the 
expression for optimal u = w + u, and solve for the trajectory 
of Z, and hence of X, but only numerically. We should also 
note that there is no guarantee in the solution above that u < 
u(t) < Vt, and < X(t) < N; this can only be verified 
numerically. 

V. Infinite-Horizon Control for Evolving Files 

We consider in this section the transmission of evolving 
files, that is files whose contents evolve and change from time 
to time. The source wishes to send the file to the destination 
and also send updates from time to time. The source need 
not know when the file changes. Updates of the file may thus 
be transmitted at times that are independent from the instants 
when the file changes. Some examples are: 

• A source has a file containing update information such 
as weather forecast or news headlines. 

• A source makes backups of some directories and store 
them at other nodes in order to improve reliability. 

The information received becomes less relevant as time 
passes. As in the original model, a relay node activates a 
time-to-live (TTL) timer when it receives a packet and deletes 
the packet when the timer expires as there is little interest in 
relaying old information. 

We are now interested in guaranteeing that there will always 
be packets in the system (as recent as possible) so that 
updated versions could be received at the destinations. The 
time horizon is now infinite so we have to restrict to those 
components of the cost defined in Section [Til] which do not 
depend on the end of the horizon. 

We wish to have X large (close to N) in order to have a 
small delivery delay (in view of (|4|). On the other hand we 
shall assign cost for low Ui to avoid old information to be 
relayed to the destination. 

We let here Z := X — N, and obtain the corresponding 
state dynamics 

dZ 

—j- = ~A out Z + w + (A in - A out )N + u 
dt 

Cost function (with Q = diag(gi, . . . , qx), <H > 0): 

POO POO 

J(oo;Z,w)= / Z T it)QZit)dt+ / w T wdt 
Jo Jo 



K rOO 

i=l J o 



z i + wf)dt 



This is a completely decoupled problem, whose solution 
involves solutions of K scalar optimal control problems. The 
i-th problem is: 



/■oo 

minJ l = / iqizf + wf)dt, 
Jo 



dzi 

= -\iz l + Wi + (fa - Aj)A^j + Ui 

where A; is the w-th element of A out , and fa is the ii-th 
element of Aj„. 

Unique stabilizing solution is (we drop the indices, and 
hence this solution is for the generic case, with everything 
being scalar): 



w 



-(pz + k), -2pX - p 2 + q = 0, 



- Afc + p(faN - XN + u) -pk = Q 

p= -A+ VA 2 +q, k=p(faN-XN + u)/y / X 2 +q 
Under this control policy, the z-th state dynamics is: 



-XiXi 



fa - pi (xi - N. t + 



+ faNi, 

+ (fa - \j)Nj 
Xi + Pi 



which is stable because Xi+pi = \J X 2 + q\ > 0. The steady- 
state value of Xi can be obtained by setting the derivative equal 
to zero in the expression above, leading to: 



(H+Pi(l- 



fa - Xi) 



:))jVi + 



Xi 



and the steady-state value of control is u°° = XiX°° — faNi . 
What remains to be shown is that there exists a choice of qi 
under which the bounds on x.i and fa are met. 

As an arbitrary special case, we picked Aj = 3, qi = 16 
(which led to pi = 2), and found that as long as fa < 3 — jf- 
both bounds are met. That is, < x°° < N% and iii < u°° < 0. 

Remark 5.1: if we want the optimal control to be linear 
(and not affine) in Xi, we start with 



Ni 



fa + (fa - Xi)Ni 



Xi +Pi 

and add to the right-hand-side the following term which is 
identically zero on the optimum trajectory at steady state 
(where on is a scalar parameter): — on{xi — xf). Now pick 
di such that all terms other than Xi are zero (and there is a 
unique such ct^), leaving us with m — ~(pi + oii)Xi. 

Remark 5.2: The above analysis can be extended to the case 
where there is a running cost on X, but then depending on the 
structure of this cost, decoupling may no longer be possible. 
Still, LQR theory would be applicable here. With coupling, 
it may not be possible to show that the bounds are satisfied 
(only through numerical computation and simulation). 

VI. Discrete-Time Control 

We shall now assume that the controlled parameters are 
updated periodically rather than continuously. Since the time 
scales involved in DTN networks are between minutes to 
hours, this is expected not to have much effect on the perfor- 
mance and yet it would decrease computing (and thus energy) 
resources. We first consider here the discrete-time version of 
the finite-horizon problem of Section llVl 



A. Finite-horizon control in discrete time 

With uniform sampling at every A units of time, the 
discrete-time version of the state equation for Z introduced 
in Section HV1 is 

Z t+ i = F e Zi + Belli + n e 

where Zt is Z(ti), with tt being the ^-th sampling time, and 
ti + i — tt — A; ue — u(t(), with control held constant over the 
subinterval [t£,tt+i)\ and Ft — $(^ + i,t^), where $ satisfies 
(and is the unique solution of) the matrix differential equation 

d$(t,T) 



dt 



A<5>{t,r) , $(t,t) = I 



that is, it is the state transition matrix associated with A. 
Furthermore, 



B e := / $(t e+1 ,T)dT B , n t := / $(^ + i,r)rfrc 

We note that if matrix A out was a constant matrix (not time 
dependent), then A would be a constant matrix, and $(t^_|_i, t£) 
would depend only on A, and likewise Bp and hi, which 
would be constants. tt) can be computed to be 



$(tt+i,tt) 
and 

B( = 
he - 



Y 







, Y := cxp (— A out A) . 



Kut[I-Y] 



rA-^AJ-A-^Z-Y] 





IA 


A-^Ai-A-yi-r]] I A 



B 



KuAi - y] 



The expression for Ft can be further simplified (approxi- 
mately) if A > is very small. To first order in A, Y = 
I — A out A, and hence 



Ft = F = 



I 
IA I 



Bt = B = 



IA 
2A-\A 



he = h = BAi n N 



Now the cost function, again as the counterpart of the one in 
Section |IV] is 

L L 

J(L; Z, u) = Zl +1 QfZ L +i+^2(ue -u) T (u e -u)+y^^ZjQZe 

i=i i=i 

where Qf is as before, and we have included an additional 
cost on intermediate values of Z, with nonnegative definite 
weighting matrix Q (which can also be taken to be zero). In 
relation to the continuous-time problem, here I is picked such 
that = t. 

As before, introducing the new control variable, wi = u^ — 
u, the state equation becomes 



Zt +1 = F e Ze + B e wt + n e , ne 



ne 



and the cost function is equivalently 

J(L;Z,u) = Zl +1 Q f Z L+1 + 



L 

wjwe 



ZjQZe 



This is a standard linear-quadratic optimal control problem, 
which admits a unique strongly time-consistent optimal solu- 
tion [14], given by 



we 



-PiSi + iF e Z? - Pe(s i+1 + S e+1 ni) 



1,2, 



where Se and s/ are generated by the recursive matrix equa- 
tions 

Pe = [I + Bj Si+iBi^Bj 



Si — Q + Fj Si + i[I — BiPeSe + i]Fe , Sl + i — Qf 



Fj[I - BiPeSi + i] T [se+i + Se+me] , sl+i 







Again, with this solution at hand, one can readily obtain the 
expression for the optimal u = w + u, and generate the 
trajectory of Z, and hence of X. 

B. Infinite-horizon control in discrete time for evolving files 

We now obtain the discrete-time counterpart of the result 
of Section [V] for the direct discrete-time version of the model 
of that section. Following the arguments of the previous 
subsection, the state dynamics in discrete time are (assuming 
that A ou t and Aj n have time-invariant, constant entries, and 5 
is as defined earlier): 

Z e+1 = GZe + Bw e + N, £=1,2,... 

where G := exp (-A out A) = diag (e~ A ' A ) 

B := A" 1 , [I - exp (-A out A)] 

N := B[(A in — A out )N + u] 

The cost function is (again Q = diag(qi, . . . , qx), Qi > 0): 



J(oo; Z, w) = (ZjQZi 



wjwe) 



K 



i=l l=\ 



This is again a completely decoupled problem, whose solu- 
tion involves solutions of K scalar optimal control problems. 
The z-th problem is: 



min Ji 
where 



-A; A 



(qiZtf + wfe), z i(e)+1 = giZu + b t w u 



X 



■0-~9i): 



bi[(jii - Xi)Ni + in 



and as before A,; is the ii-th element of A out , and fa is the 
M-th element of A jn . 

The unique stabilizing solution is (we drop the indices, and 
hence this solution is for the generic case, with everything 
being scalar): 

b 



w = -PSgz - P(s + Sn), P = 



l + Sb 2 ' 

S = q + g 2 S[l- bPS], s = g[l- bPS] [s - Sn], 
1 



Sn 



l + b 2 S 



-Si 



l + b 2 S- g 

where it can be shown that 1 + b 2 S — g ^ 0, and hence the 
expressions for S and s are well defined. The control policy 
simplifies to 

bSg bS 

w ~ ^ITsp z ^ T~g~T¥s n 

which is again for the generic case (for the i-th control all 
quantities above will have the index i). Under this control 
policy, the i-th state dynamics is: 

Qi 1 



Z i(t+1) 



1 + Sib 2 



1 



SA 2 ' 



which is stable because < gi/(l + Sib 2 ) < 1, with the reason 
being that < gi < 1 and 1 + Sibf > 1. The steady-state value 
of Zi is: 



(1 



where n, 



A, 



(l-e x > A ){(fa-\ l )N + ' 



and the actual state is of course irj 



N. The steady- 



state open-loop value of the optimal control is 

hSi 



hSm 



■n, 



If A is very small, asymptotically Si A 
and the optimal control policy becomes 



opt 

U; = — 



A,, _ 



(1 



V qi + K 



f )(fa - \) + \i - yt+A? 



N 



Again in the asymptotic case (for small A), the steady-state 
value of the state (i-th component) is 



N ■ 



Xi(fU-Xi)N 
Note that as A — > 0, x°° — > N. 



A- 



A, 



y/n + x i 



zUiA 



S = 



2b 2 



b 2 q + g 2 - 1 + ^J{b 2 q + g 2 - l) 2 + 4qb 2 



VII. Numerical Results 

Here we present numerical results for the dynamic control of 
timers; in particular, we report on the dynamics of the system 
in some reference cases. The numerical integration procedure 
was performed as follows 

1) the solution P of the differential matrix Riccati equation 
is calculated over the interval [0, t/]; 

2) the dynamics of k have been obtained integrating back- 
wards the related differential equation; 

3) finally, an interpolated version of vector k and matrix P 
are input to the forward integration of the ODE describing 
the state trajectory X; 

4) the control law u and the corresponding message delivery 
probability are derived accordingly. 
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Fig. 1: The minimum value C4/C3 as a function of ci/cj, such 
that R > 0; the solid line refers to the uniform case = 
(1, 1, l)/\/3, A out = diag{(l, 1, 1)}, the dot-dashed line to the 
non-uniform case = (1, 3, 5)/\/35, A out = diag{(l, 2, 4)}. 



For the cases described in the following we employed the 
Matlab ODE suite. Mobile nodes intermeeting intensities are 
dimensioned using the Random Waypoint (RWP) synthetic 
mobility model for which intermeeting intensities can be 
calculated [9]. The reference intermeeting intensity, for the 
following experiments, is Ao = 2.7875 ■ 10~ 4 ; this value 
is experienced by RWP mobile nodes moving on a square 
playground of side L = 1000 m, speed v — 4 m/s and radio 
range 20 m. 

Choice of parameters of the linear quadratic optimization 

We first study the parameters of the optimization problem 
and related constraints. In particular, Ci,c 3 and C4 enable 
to tune the optimization based on specific weighs that we 
can assign to the delay, energy and the memory; we notice 
that ci , C3 and C4 appear in the definition of R, whereas ci 

■ , tt ■ • rr ( A -BB T /C 2 ' 

appears in the Hamiltonian matrix ti = I ^ ^j, 

Note also that we had set C2 = 1 at the outset, as part of 
normalization of the cost. 

In particular, c\ , C3 and C4 are such that 

R = -ciK d K* d + c 3 I + Ci A 2 out > 

In view of the expression for R, for a given value of C1/C3, 
a minimum value for C4/C3 exists such that R > 0. Here, we 
provide a simple sufficient condition. 
Proposition 7.1: Let 



K 



A 4 

1 ^id 



if C4 > a(ci — c 3 ), then R > 0. 

The proof of the above statement is reported in the Appendix. 
We observe that for c\ < C3 and C4 > the condition 
is satisfied. We would also like to give some insight in the 
relative dimensioning of C1/C3 and C4/C3 for the outer region 
ci > C3. We derived numerically the minimum values of C4/C3 
above for which R > in case K = 3. As depicted in Fig. [T] 



we reported both on the case of uniform entries for A^ and 
A out and the case when they are different; we referred to the 
normalized case when ||Ad|| = 1 for the sake of the clearness. 
There, we can distinguish the inner region for c%/cs < 1 
where Prop. l7.ll holds. and the outer region where C4/C3 must 
increase in order to ensure R > 0; for the tested values the 
graph shows a quasi linear shape and the lower bound proves 
tight in the case of uniform entries (L.B. lines in Fig. [T). 

This means that at the increase of C1/C3, i.e., for larger 
relative weight given to the delivery probability in the cost 
function, C4/C3 has to increase: the relative weight given to 
the energy cost has to increase in order to maintain the problem 
positive definite: this corresponds to the intuition that it is not 
possible to overweight the delivery probability of a message 
against the energy expenditure. 

The effect of the energy constraint ( C4) 

First, we describe the impact of the constraint on energy 
C4; in particular, we considered terminal time r = 3600 s. 
Ni - 1 = N 2 = N 3 = 50. Also, we let u = 0. 

The coefficients of the optimization problem are c\ = 1/Aq 
and C3 = 1. We considered uniform intermeeting intensities, 
A si = Xjd = An, for i = 1,2, 3. 

We compared the performance of the system for three values 
of C4: 0.05, 0.01 and 0.005. As expected, see Fig. [2] when 
c 4 = 0.005, i.e., the constraint on the energy expenditure for 
message transmission is smaller, the effect of timers is milder. 
For this setting, the message delay CDF reaches the unitary 
value much before the terminal time. The controlled dynamics 
of the number of infected nodes are monotonic as in the case 
of the plain, uncontrolled, two-hop routing. 

Conversely, for larger values, i.e., C4 = 0.01, the message 
delay CDF reaches the unitary value slightly later than in 
the previous case (see Fig. |2^), since the effect of timers 
becomes more relevant (see Fig. Furthermore, the change 
of convexity of the infected nodes dynamics (see Fig. Hh) is 
marked. For C4 = 0.05, the effect of the larger weight given 
to the energy results in a non-monotonic dynamics in the 
number of copies in the system. We notice that, at the opposite 
ends, the stronger constraint on the energy leads to a much 
smaller number of copies in the system, 10 per class in case 
of C4 = 0.005, and 4 in case of C4 = 0.05: nevertheless in the 
latter still D(t) ~ 1. This first result already indicates that via 
proper parameter weighting we can achieve consistent savings 
of network resources. 

The case of inhomogeneous classes 

So far we investigated the properties of the system in the 
case when the classes of mobiles are homogeneous. Here, 
Aid = A2d/2 = 3A3d/2, whereas the intermeeting intensity 
with the source A S i = Ao, i=l,2,3. In order to maintain consis- 
tency with the previous cases, we normalized the intermeeting 
intensities ||A OMt || = ||Aj n ||. Also, in order to meet the 
constraints, for this choice of the parameters, a suitable setting 
is c 4 = 0.0025, c 3 = 1/2, c x = 1/2A§, and = -0.004, 
i = 1,2,3. 





In this numerical evaluation, see Fig. [4] the control dis- 
criminates classes with higher chances to deliver the message 
(classes 2 and 3) and the first class; hence, higher timer 
rates are assigned to the first class to limit the number of 
copies forwarded to nodes with lower intermeeting rates to 
the destination. 

In addition to the setting of Fig. |4] we considered also 
different intermeeting intensities with the source at different 
classes (see Fig \5}. In addition to the setting described in the 
previous case, in particular, A s i = 0.7A S 2 and A S 3 = 1.3A S 2. 
As depicted there, the change of the relative intermeeting rates 
within the three classes has a marked impact into the way the 
control is performed compared to the case of Fig. |4] In this 
case, in fact, the timeout rate, with respect to classes 2 and 3 
is still larger, in order to limit the increase of the number of 
messages; notice, though, that the infected nodes dynamics of 
those two classes is basically the same of the previous case. 
But, the control of the timeouts for class 1 is much milder than 
seen previously, due to the lower intermeeting rate within that 
class, which reduces the need for high timeout rates. 



A. The use of reference timeouts rates 

In the last two cases, i.e., Fig. ®> and Fig. \5]p, we expanded 
the time scale around r in order to confirm the numerical 
stability of our solution. In particular, the final value of 
the control is dictated by the reference value u. This also 
suggests that finer tuning of the control can be obtained using 
different values for each Ui. With respect to Fig. |6j?, we 
observe that this fine tuning is beneficial: under the same 
settings of Fig. [5J5, the system is able to reach the desired 
high delivery probability, whereas the number of copies in the 
system decreases compared to the use of uniform penalty on 
timeout rates. 

Remark 7.1: The numerical evaluations represented before 
do not exhaust the range of the possible parameters of the 
problem. For the choices showed before, we limited to cases 
where the optimal solution is compatible with the constraints 
< Xi < N and Ui < 0; for cases of practical 
interest, due to the particular structure of the two-hop routing 
protocol, the upper bound on the number of nodes is usually 
satisfied. We found that a crucial parameter is the reference 




t [s] t [s] 

(a) (b) 

Fig. 4: Dynamical control of the timers, non-uniform case, N — 151, r = 3600 s; (a) delay CDF (upper) and dynamics 
of infected nodes per class Xi(t) (lower); (b) control Ui. Parameters are C4 = 0.0025, c 3 = 1/2, c 2 = 1 and ci = l/2Ag. 
Intermeeting intensities with destination A^ = l/2A2d = i/2\zd, ||A ou t|| = ||Am|| = Aq. 




Fig. 5: Dynamical control of the timers, non-uniform case, N — 151, r = 3600 s; (a) delay CDF (upper) and dynamics of 
infected nodes per class Xi(t) (lower); (b) control Uj. Parameters are as in Fig. [4] A s i = 0.7A S 2 and A S 3 = 1.3A S 2; again 
||A out || = ||A m || = A . 
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Fig. 6: Dynamical control of the timers, non-uniform case, N — 151, r = 3600 s; (a) delay CDF (upper) and dynamics of 
infected nodes per class Xi(t) (lower); (b) control Uj. Parameters are as in Fig. [5] u = — (4.1, 5.1, 6) T x 10 -3 . 



intermeeting intensity Ao, which appears to strongly impact 
the sensitivity of the optimal solution; in particular, it may 
determine settings where the constraint Xi > cannot be 
attained. More precisely, the approach seems to show some 
limitation when dealing with small values of Ao- 

VIII. Conclusion 

We have focused in this paper on controlling the spreading 
of message under two-hop relay forwarding. We exploited the 
linear form of the dynamics of this regime to study the control 
problem within the linear quadratic control framework. This 
allowed us to study numerically the tradeoff achievable by 
tuning various parameters that define the cost. 
There exist several aspects that were not covered in the present 
work, which deserve further investigation. We showed that 
there exist constraints on the choice of optimization param- 
eters; this relates to the need to keep the cost function well 
defined and bounded. We also remarked the need to satisfy the 
constraints on the dynamics of the number of infected nodes: 
to this respect, the numerical evaluations reported before do 
not exhaust the range of the possible parameters. Due to 
the particular structure of the two-hop routing protocol, the 
upper bound on the number of nodes is usually satisfied. 
We observed, though, that the reference intermeeting intensity 
Ao strongly affects the sensitivity of the optimal solution; 
in particular, it may determine settings where the constraint 
Xi > cannot be attained, especially for very small values 
of Ao. Notice, however, that even when the constraints are 
satisfied for the dynamics of the averages, they need not be 
satisfied each sample. An interesting question is to determine 
the probability that they are not satisfied by some realization. 
From a practical standpoint, there are other directions for 
future work. The implementation of the control in two-hop 
routing can be done at the source only, and timers can be reg- 
ulated through appropriate time-stamping of message copies 
(relays perform message discarding simply comparing the time 
elapsed since message generation and the control reported on 
the message copy). The crucial problem, conversely, is the 
estimation of the system parameters, namely, N and \j at 
the source node. Also, previous work [20], showed that in 
the 1-dimensional case, when these parameters are unknown, 
it is still possible to obtain a policy that converges to the 
optimal one by using some auto-tuning mechanism. Using 
stochastic approximation this policy is shown to be optimal 
for the average cost criterion. 

An interesting research direction would be to combine the 
optimal control techniques presented in this paper with coding 
techniques such as fountain codes or network coding [21], 
[22]. 
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Appendix 

Deriving the success probabilities 

We follow the derivation in [16, Appendix A] which we 
extend so as to handle the multi-class case, as well as to handle 



the case of non-homogeneous parameters (we allow 

Pr(T d >t + h\T d >t,3 t )=l-hY J fc(*)A«(t) + o(h) 

i 

which implies 

Pr(T d > t + h\J t ) = Pr(T d > t|J t )(l - ft^&W-M*)) 
As = 1 - Pr(T d > t\5 t ) we get 



Thus 



d£(j 



1 - *(t) 



= ^&(t)A id (t) 



log(l - (t)) = 2 A M / + c i 



and hence 



= 1 - C 2 exp ( - V A i4 / £ (s)d 



(8) 



For < large, ^ (t) should tend to one, which it does. Also as 
t — ► 0, tf(i) should go to zero. This last condition implies that 
C 2 = 1. 



Proof of Prop. \7.1\ 

Proof: R > if and only if -ci/c 3 A d Aj + 7 + 
C4/c3A^ nf > 0. Also, ||Arf|| 2 is the only positive non-zero 
eigenvalue of M = A d A d , whose eigenvector is Ad. M is 
symmetric, so let V an orthogonal matrix such that 

diag(||A d || 2 ,0,...,0) T = F T M^, 

where in particular V — (Aj, of, ■ • • , c^) and Cj € ker(Af). 
Also, R > iff R = V T RV > 0, i.e. 

R = I — -diagfl | A d | | 2 , 0, . . . , 0) + -V T A out A out 7 > 0. 

C3 C 3 

Let ei = (1,0,.. .,0): the sufficient condition is obtained 
since if e[Rei > 0, then vRv > for Vu £ R fe . Hence, 

dReJ = (l-^)||Ad|| 2 + ^||A out y ei || 2 
= (i-£l)||Ad|| 2 + ^||A OTf Ad|| 2 

from which the statement follows. ■ 

Calculation of P 

In the following we derive a closed form solution for P; we 
assume that (A ou t)u > 0, for i = 1,2,..., aH The solution 
of the matrix Riccati equation is given by P = P2P1 1 , where 
Pi , P2 are 2K x 2 A matrices solutions of 



Pi(t f ) = I 2K , P 2 (t f ) = Q f , (9) 



2 The expression that we derive in the following requires A «t to be 
invertible but not diagonal. 



where I 2K is the 2A x 2A identity matrix; Pi is guaranteed 
to be invertible for all < t < tf. In our case, if we choose 
C2 = 1, the Hamiltonian matrix is 



H 



A -BB T 
-A T 



V 



A ou t 





-I K 




I K 

















A ou t 


-Ik 











J 



The associated dynamical system (O solves for 



Pi(t) 



= e H(t-t f : 



Q f 



so that we are interested in the explicit calculation of the 

matrix e H ( t ~ t - f ' ) — YlkLo ^~kP H k where for the sake of 

notation, x = (t — tf). 

The fc-th power of H can be derived as follows 
Proposition 9.1: H° = I AK , H 1 = H, for k > 1: 



H k = 



H K = 



-A k 
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1Y out 
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if fc is odd 



if k is even 



The above formula can be easily verified by induction from 
( fTOl i; hence, if By is the «j-th K x K block of e Hx , i,j = 
1, 2, 3, 4, we obtain 



Ex 
fe! 

k=Q 



where ( J ff fc )y the K x K ij-block of H k , i, j = 1, 2, 3, 4. 

From (ED, £?i2 = £31 = £32 = £41 = E 42 = E A3 = 0; 
the diagonal entries are En = E^ 3 = X)^Lo ' ^~ A " ' ' ~~ 



e A °" ta: whereas E22 = E44 = I K . Non-zero off-diagonal 
entries require some calculations, which bring 



A'! 



£21 



<ut X 1 



^-out ( I'- 



34 



-A; 



-1 (pKutX _ j \ 

•ut\ C 1 Kj- 



For ease of presentation, given diagonal matrix A, we define 
6(A) and 8(A) the diagonal matrix such that G(A)u — 
cosh(A,i) and §(A)u = sinh(A,i), respectively. Hence, it 
follows 



E 



13 



-A * t S(A 0Ut x), E 2i = A otl 3 t (S(A ouf x) - xA out ) 



E u = A - u 2 (e(A out x)-I K ), E 23 = -A-* t (e(^outx)-I K ) 
Thus, we obtain 



P2 = 
Pi = 



-A: 



R 



R 



-A ou tx 



-A out x 



A-^(e(A out x)-I K )R 
J K + Kut(H A outx) - xA out )R 



The inverse of Pi can be obtained leveraging the block form 
and requiring Pi-Pf 1 = I 2K 



A 3 = M-\x)A^ t (I K - e A — ), A 4 = M~\x) 

where M x = I K + A m ; 3 t [(§(A out a;) - xK out )-{Q{K out x) - 
I K ){e K ^ x -I K )]R. 



Pu = A-^-e^nM-'R 
P 21 = i?M- 1 A OM 1 4 (/ K - e A — ) 
P 22 = RM- 1 



We notice that matrix M x is symmetric and invertible; in 
fact, let f(x) — sinh(x) — x+(cosh(x) — 1)(1 — e x ); it follows 
that / = ~(e x -l) 2 , so that [(§(A out x)-xA out )-(e(A out x)- 
I K )(e A ° utX -I K )] > since -t f <x<0. Then, v T M (x)v > 
for Vu 7^ 0, so that M{x) > 0. P is symmetric, as expected. 




from which we obtain 



Finally, we obtain 



